Spatial data in R: an introduction to the sf package

---

Installation and loading

To install sf from CRAN:

install.packages("sf")  

or the development version from GitHub:

devtools::install_github("edzer/sfr")

then load:

library(sf)

NB MacOSX and LINUX users need to install a number of geospatial libraries (GEOS, GDAL, and proj.4).

The tidyverse package also needs to be loaded:

library(tidyverse)

Example data

A GeoJSON of Greater Manchester’s wards was created from a vector boundary file available from ONS’s Open Geography Portal. The GeoJSON is projected in British National Grid (EPSG:27700) and originally derives from the Ordnance Survey.

Point data are supplied by data.police.uk and represent incidents of anti-social behaviour and crime recorded by the Greater Manchester Police during February 2017. The incidents are supplied with latitude and longitude coordinates which have undergone an anonymisation process.

Reading and writing spatial data

Reading spatial data (polygons)

bdy <- st_read("data/wards.geojson", quiet = TRUE)

Reading data with coordinates (points)

pts <- read_csv("data/2017-03-greater-manchester-street.csv")
pts <- st_as_sf(pts, coords = c("Longitude", "Latitude"))

Writing spatial data

st_write(bdy, "boundaries.shp")

sf objects

Attribute table (dataframe) AND geometry (list-column) with coordinates, CRS, bbox

class(bdy)

[1] "sf"         "data.frame"

head(bdy)

Simple feature collection with 6 features and 12 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 351664 ymin: 394521.5 xmax: 369078.9 ymax: 407481
epsg (SRID):    27700
proj4string:    +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs
  objectid    wd16cd      wd16nm wd16nmw   lad16cd lad16nm  bng_e  bng_n     long      lat
1      772 E05000851        Ince    <NA> E08000010   Wigan 359959 405278 -2.60570 53.54262
2      773 E05000852  Leigh East    <NA> E08000010   Wigan 367298 400695 -2.49447 53.50194
3      774 E05000853 Leigh South    <NA> E08000010   Wigan 366260 398660 -2.50990 53.48358
4      775 E05000854  Leigh West    <NA> E08000010   Wigan 363428 400762 -2.55282 53.50228
5      776 E05000855 Lowton East    <NA> E08000010   Wigan 362642 397699 -2.56431 53.47470
6      777 E05000856      Orrell    <NA> E08000010   Wigan 353321 404083 -2.70568 53.53133
  st_areasha st_lengths                       geometry
1    4780173   11442.63 POLYGON((359198.460435328 4...
2    4972092   11056.20 POLYGON((367700.613520717 4...
3    6521214   12377.97 POLYGON((364489.980933077 3...
4    7315344   16159.92 POLYGON((363107.910792193 3...
5   11185415   17990.73 POLYGON((361792.312268886 3...
6    8884317   19650.85 POLYGON((353418.76434581 40...

as.tibble(bdy)

bdy_df <- st_set_geometry(bdy, NULL) # remove geometry
head(bdy_df)

st_geometry(bdy) # print geometry

Geometry set for 215 features 
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 351664 ymin: 381168.6 xmax: 406087.5 ymax: 421039.8
epsg (SRID):    27700
proj4string:    +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs
First 5 geometries:

Manipulate using dplyr

Rename variables

bdy <- bdy %>% 
  select(ward = wd16nm, census_code = wd16cd, borough = lad16nm) # rename variables
glimpse(bdy)

Observations: 215
Variables: 4
$ ward        <fctr> Ince, Leigh East, Leigh South, Leigh West, Lowton East, Orrell, Pemberto...
$ census_code <fctr> E05000851, E05000852, E05000853, E05000854, E05000855, E05000856, E05000...
$ borough     <fctr> Wigan, Wigan, Wigan, Wigan, Wigan, Wigan, Wigan, Wigan, Wigan, Wigan, Wi...
$ geometry    <simple_feature> POLYGON((359198.460435328 4..., POLYGON((367700.613520717 4......

Count frequency of wards by borough

bdy %>% 
  group_by(borough) %>% 
  count() %>% 
  arrange(desc(n)) %>% # sort in descending order
  st_set_geometry(., NULL) # hide geometry

Select features

chorlton <- bdy %>% 
  filter(ward == "Chorlton")
plot(st_geometry(bdy))
plot(st_geometry(chorlton), col = "red", add = TRUE)

Using sf functions to add geometry column to dplyr chain

bdy <- bdy %>% 
  mutate(area = st_area(.)) # returns the area of a feature
bdy %>% 
  select(ward, area) %>% 
  arrange(desc(area)) %>% 
  slice(1:10) %>% 
  st_set_geometry(., NULL)

Projection

Check and assign CRS

st_crs(pts) 

$epsg
[1] NA

$proj4string
[1] NA

attr(,"class")
[1] "crs"

pts <- st_set_crs(pts, 4326) # assign Lat/Long (epsg:4326)
st_crs(pts)

$epsg
[1] 4326

$proj4string
[1] "+proj=longlat +datum=WGS84 +no_defs"

attr(,"class")
[1] "crs"

Reproject CRS

bdy_WGS84 <- st_transform(bdy, 4326)
st_crs(bdy_WGS84)

$epsg
[1] 4326

$proj4string
[1] "+proj=longlat +datum=WGS84 +no_defs"

attr(,"class")
[1] "crs"

Convert to and from sp objects

Convert to and from sp objects

bdy_sp <- as(bdy, 'Spatial')
class(bdy_sp)

[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

bdy_sf <- st_as_sf(bdy_sp)
class(bdy_sf)

[1] "sf"         "data.frame"

Spatial operations

Buffer features

buffer <- chorlton %>% 
  st_buffer(dist = 1000)
plot(st_geometry(buffer))
plot(st_geometry(chorlton), col = "red", add = TRUE)

Buffer and intersect

pts_sub <- bdy %>%
  filter(ward == "Chorlton") %>%
  st_buffer(dist = 1000) %>%
  st_intersection(st_transform(pts, 27700)) # reproject pts to BNG (epsg:27700)
plot(st_geometry(buffer))
plot(st_geometry(chorlton), col = "red", add = TRUE)

plot(st_geometry(pts_sub), col = "black", add = TRUE)

pts_sub %>% 
  group_by(`Crime type`) %>%
  count() %>% 
  arrange(desc(n)) %>% 
  st_set_geometry(., NULL)

Points in polygon

pts %>% 
  filter(`Crime type` == "Vehicle crime") %>%  
  st_join(bdy_WGS84, ., left = FALSE) %>% 
  count(ward) %>% 
  arrange(desc(n)) %>% 
  st_set_geometry(., NULL)

bdy_pts <- pts %>% 
  filter(`Crime type` == "Vehicle crime") %>%  
  st_join(bdy_WGS84, ., left = FALSE) %>% 
  count(ward)

Plotting

Base plots

plot(bdy_pts) # plots small multiples if dataframe has several attributes

plot(bdy_pts["n"]) # select the appropriate attribute to plot a single map

library(RColorBrewer) ; library(classInt)
pal <- brewer.pal(5, "RdPu")
classes <- classIntervals(bdy_pts$n, n=5, style="pretty")$brks
plot(bdy_pts["n"], 
     col = pal[findInterval(bdy_pts$n, classes, all.inside=TRUE)], 
     main = "Vehicle crime in Greater Manchester\nMarch 2017", axes = F)
legend("bottomright", legend = paste("<", round(classes[-1])), fill = pal, cex = 0.7) 

ggplot2

devtools::install_github("tidyverse/ggplot2") # NB need development version for geom_sf()

ggplot(bdy_pts) +
  geom_sf(aes(fill = n)) +
  scale_fill_gradientn('Frequency', colours=RColorBrewer::brewer.pal(5,"RdPu"), 
                       breaks = scales::pretty_breaks(n = 5)) +
  labs(fill = "Frequency",
       title = "Vehicle crime",
       subtitle = "March 2017",
       caption = "Contains OS data © Crown copyright and database right (2017)") +
  theme_void()

plotly

library(plotly)
p <- ggplot(bdy_pts) +
  geom_sf(aes(fill = n, text = paste0(ward, "\n", "Crimes: ", n))) +
  scale_fill_gradientn('Frequency', colours=RColorBrewer::brewer.pal(5,"RdPu"), 
                       breaks = scales::pretty_breaks(n = 5)) +
  labs(fill = "Frequency",
       title = "Vehicle crime",
       subtitle = "March 2017",
       caption = "Contains OS data © Crown copyright and database right (2017)") +
  theme_void() + 
  coord_fixed(1.3)
ggplotly(p, tooltip = "text")

leaflet

library(leaflet)
pal <- colorBin("RdPu", domain = bdy_pts$n, bins = 5, pretty = TRUE)
leaflet(bdy_pts) %>% 
  addTiles(urlTemplate = "http://{s}.tiles.wmflabs.org/bw-mapnik/{z}/{x}/{y}.png",
    attribution = '&copy; <a href="http://www.openstreetmap.org/copyright">OpenStreetMap</a>, <a href="https://www.ons.gov.uk/methodology/geography/licences">Contains OS data © Crown copyright and database right (2017)</a>') %>% 
  addPolygons(fillColor = ~pal(n), fillOpacity = 0.8,
              weight = 1, opacity = 1, color = "black",
              label = ~as.character(ward)) %>% 
  addLegend(pal = pal, values = ~n, opacity = 0.7, 
            title = 'Vehicle crime (March 2017)', position = "bottomleft") %>%
  addMiniMap(tiles = providers$CartoDB.Positron, toggleDisplay = TRUE)